function [ outdata ] = ratheartdrawallspio( basedir, subdirs, overlayrange, varargin )
% ratheartsegmentation: draw the spio and a myocardium region around it; 
% draw graphs over time
%
% Myocardium segmentation of the left ventricle for short axis MR scans.
% Generates a curve showing the evolution of R2 values within each region
% over time.
%
% Usage: 
%    outdata = ratheartdrawspio(basedir, subdirs);
%
% where:
%    basedir = the dicom directory of a specific rat
%    subdirs = a cell array of strings with subdirectories pointing to the
%              T2 subdirectories.
%    outdata = a struct array, containing for each day some relevant data.
%              outdata(day).R2map contains the registered R2map for
%              example.

% Adapted by Henk Smit, based on adaptation by Stefan Klein of code by Gabrielle Giorico, September 2011.


%force redrawing ROIs
force = false;
if nargin > 3
    force = varargin{1};
end

% determine postproc dirs.
nb_days = numel(subdirs);
for day = 1:nb_days
  studydir = fullfile(basedir, subdirs{day} );
  % remove slashes at the end
  while ( studydir(end) == '/' || studydir(end) == '\' )
    studydir = studydir(1:(end-1));
  end
  outdir{day} = [ studydir 'postproc\' ];
  if ~exist(outdir{day}, 'dir')
    error( 'ratheartsegmentation:invalidoutputdir', '%s : invalid postproc dir. Run first ratheartT2fit!', outdir{day});  
  end
end

% all results are written to the first day's postproc directory
resultdir = outdir{1};

%Original image
firstTEinfo = dicominfo(fullfile(resultdir,'FirstIm.dcm'));
firstTE=double(dicomread(firstTEinfo));

% Load the first day's R2 map
info_zeroline = dicominfo( fullfile( resultdir, 'R2map.dcm' ) );
R2mapfirstday = double(dicomread(info_zeroline));


% Initialization
array_meanr2 = []; 
array_std = [];
array_crlb = [];

    % Loop over all the days
for day=1:nb_days
    
    %Load maps and images

    infoR2map = dicominfo( fullfile( outdir{day}, 'R2map.dcm') );
    infoR2errormap = dicominfo( fullfile( outdir{day}, 'R2errormap.dcm') );
    infoTE = dicominfo( fullfile( outdir{day}, 'FirstIm.dcm') );

    R2map = double(dicomread(infoR2map));
    R2maps(day).image = R2map;
    R2errormap= double(dicomread(infoR2errormap));
    TE = double(dicomread(infoTE));
    MRIs(day).image = TE;
    mask = int16(R2maps(day).image~=0);
    
    %when ROI is to be drawn here:

        % size of figure for manual segmentation:
        figuresize = [100 100 1000 800];
        spioroifile = [outdir{day} 'spioroi.mat'];
        myocardroifile = [outdir{day} 'myocardroi.mat'];
    if ( ~exist( spioroifile, 'file' ) ||  ~exist( myocardroifile, 'file' ) || force )

        figure( 'position', figuresize );  
        imshow(MRIs(day).image,[]);    
        title( 'draw a region with SPIO (end with double-click)' );
        drawnow;
        SpioMask = roipoly;
        title( 'draw a region with myocardium (end with double-click)' );
        M=roipoly; % outer and inner circle of myocardium
        %M = M - N;
        close;
        MyocardMask = M-SpioMask;
        map = uint16(R2maps(day).image);
        R2mapmask=logical(map./map);
        MyocardMask = MyocardMask.*R2mapmask;
        save( spioroifile, 'SpioMask');
        save( myocardroifile, 'MyocardMask');
    else
        if day==1
        disp( 'using previously selected SPIO and Myocardium. Use "ratheartdrawspio ( basedir, subdirs, true )" in ratheartpostproc.m to redraw the ROIs.' )      
        end
        load( spioroifile, 'SpioMask');
        load( myocardroifile, 'MyocardMask');          
    end

    masks(day,1).image = SpioMask;
    masks(day,2).image = MyocardMask;
    spioedges(day).image=edge(uint16(SpioMask),'canny');


    % To find the values corresponding to a region in the R2map
    vec_region_day = [];
    vec_std_day = [];
    vec_crlb_day = [];

    for r=1:2  %spio and myocardium
        r2values = R2map( (masks(day,r).image == 1) & (R2map>0) ) ;
        r2values_final = r2values/10000;
        r2errvalues = R2errormap( (masks(day,r).image == 1) & (R2map>0) );
        r2errvalues_final = r2errvalues/10000;

        % store mean r2, std r2, and mean crlb
        vec_region_day = [vec_region_day, mean(r2values_final)];
        vec_std_day = [vec_std_day, std(r2values_final)];
        vec_crlb_day = [vec_crlb_day, median(r2errvalues_final)]; 
       
    end % end of the region loop
    
    % Update the matrix containing the mean value of the region per day
    % (each line is a day)
    array_meanr2 = [array_meanr2; vec_region_day];
    array_std = [array_std; vec_std_day];
    array_crlb = [array_crlb; vec_crlb_day];
    
    % store for user:
    outdata(day).R2map = R2map;
    outdata(day).R2errormap = R2errormap;
    outdata(day).firstTE = TE;
    outdata(day).postprocdir = outdir{day};
    outdata(day).segmentation = mask;
    
end % End of the day loop

% Plot
figure('Name', 'Analysis results', 'position', [50 400 800 500]);
vec_xaxis = [];

% TO DO: Look for the command lign that create the vector directly
for i=1:nb_days
    vec_xaxis = [ vec_xaxis, i];
end

% Graphs
shiftbetweenstdandcrlb = 0.05;
ylim = [0 0.15];
labels = {'SPIO','Myocardium'};
for r=1:2
    subplot(1,2,r);
    errorbar( vec_xaxis, array_meanr2(:,r)', array_std(:,r), '.b-'); 
    hold on;
    if r==1
        % dummy, to fake legend
        errorbar( 10, 1, 1, '.r-'); 
    end
    errorbar( vec_xaxis + shiftbetweenstdandcrlb, array_meanr2(:,r)', array_crlb(:,r), '.r'); 
    axis([ 0.8 nb_days+0.2 ylim]);
    h = title(labels(r));
    pos = get(h,'Position');
    pos(2) = ylim(1)+ (ylim(2)-ylim(1))*0.85;
    set( h, 'Position', pos);
    xlabel('Scans');
    set(gca,'XTick', 1:nb_days );
    ylabel('R2 values [ms]');
    if r==1        
        lh=legend({'st.dev.', 'CRLB'}, 'Location', 'North','Orientation', 'horizontal');
    end
end

pos = get(lh, 'Position');
pos([1 2]) = [0.4 0];
set(lh, 'Position', pos);

%load parameters if exist
if exist(fullfile(basedir,'Estimated_parameters.mat'))==2
    load(fullfile(basedir,'Estimated_parameters.mat'))
end 



%store and write all results
outputfile=[resultdir 'results.txt'];
fid = fopen(outputfile, 'a');
% fprintf(fid, 'T2_Spio T2CR_Spio T2STD_Spio T2_Myocard T2CR_Myocard T2STD_Myocard R2_Spio R2CR_Spio R2STD_Spio R2_Myocard R2CR_Myocard R2STD_Myocard\n');
fprintf(fid, 'Day T2_Spio T2_Myocard R2_Spio R2_Myocard R2CR_Spio R2STD_Spio  R2CR_Myocard R2STD_Myocard\n');
for day=1:nb_days
    fprintf(fid, '%i %f %f %f %f %f %f %f %f\n', day,1/(array_meanr2(day,1)),1/(array_meanr2(day,2)),array_meanr2(day,1),array_meanr2(day,2),array_crlb(day,1),array_std(day,1),array_crlb(day,2),array_std(day,2) );
    figure('Name', ['Overlay Day ' num2str(day)], 'position', [800+(day*150) 600-(day*150) 500 500]);
    sc(R2maps(day).image,overlayrange, jet,sc(MRIs(day).image, [0 2000],gray),(masks(day,2).image+3.*masks(day,1).image)<1,sc(spioedges(day).image,[1 1],gray),spioedges(day).image==1);
    colorbar('south','XColor','white');
    parameters(1,day,1) = array_meanr2(day,1);
    parameters(1,day,2) = array_meanr2(day,2);
    parameters(2,day,1) = array_crlb(day,1);
    parameters(2,day,2) = array_crlb(day,2);
    parameters(3,day,1) = array_std(day,1);
    parameters(3,day,2) = array_std(day,2);

end

    save([basedir, filesep,'Estimated_parameters.mat'],'parameters')
    disp( ['Saved results to ' outdir{1} 'results.txt'] );
    fclose(fid);










   